"""
CBMPy: Constraint Based Modelling in Python (http://pysces.sourceforge.net/cbm)
============
Copyright (C) 2009-2012 Brett G. Olivier, VU University Amsterdam, Amsterdam, The Netherlands

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>

Author: Brett G. Olivier
Contact email: bgoli@users.sourceforge.net
Last edit: $Author: bgoli $ ($Id: CBCPLEX.py 990 2012-02-15 13:52:11Z bgoli $)

"""
import os, time, cStringIO, gc, random
import numpy
import cplex
import CBWrite, CBTools
from CBConfig import __CBCONFIG__ as __CBCONFIG__
__DEBUG__ = __CBCONFIG__['DEBUG']
__version__ = __CBCONFIG__['VERSION']

##  class FluxSensitivity(object):
    ##  """
    ##  Contains variable (Flux) sensitivities and reduced cost
    ##  """
    ##  flux = None    
    ##  lower_bound_low = None
    ##  lower_bound_high = None
    ##  lower_bound_value = None
    ##  upper_bound_low = None
    ##  upper_bound_high = None
    ##  upper_bound_value = None
    ##  objective_coefficient_value = None
    ##  objective_coefficient_low = None
    ##  objective_coefficient_high = None
    ##  reduced_cost = None
    
    ##  def getUpperBoundRange(self):
        ##  """
        ##  Returns the absolute range of the upper bound.
        ##  """
        ##  return abs(self.upper_bound_high-self.upper_bound_low)

    ##  def getLowerBoundRange(self):
        ##  """
        ##  Returns the absolute range of the lower bound.
        ##  """
        ##  return abs(self.lower_bound_high-self.lower_bound_low)

    ##  def getObjectiveCoefficientRange(self):
        ##  """
        ##  Returns the absolute range of the objective flux coefficient.
        ##  """
        ##  return abs(self.objective_coefficient_high-objective_coefficient_low)

##  class RHSsensitivity(object):
    ##  """
    ##  Class containing RHS sensitivities
    ##  """
    ##  species = None
    ##  rhs_low = None
    ##  rhs_high = None
    ##  rhs_value = None
    
    ##  def getRHSRange(self):
        ##  """
        ##  Returns the absolute range of the objective flux coefficient.
        ##  """
        ##  return abs(self.rhs_high-rhs_low)

##  class Sensitivities(object):
    ##  """
    ##  Container class that holds sensitivities
    ##  """
    ##  flux_ids = None
    ##  species_ids = None
    
    ##  def __init__(self):
        ##  self.flux_ids = []
        ##  self.species_ids = []
        
    ##  def addFluxSensitivity(self, fsense):
        ##  self.__setattr__(fsense.flux, fsense)
        ##  self.flux_ids.append(fsense.flux)
        
    ##  def addRHSsensitivity(self, rhsense):
        ##  self.__setattr__(rhsense.species, rhsense)
        ##  self.species_ids.append(rhsense.species)
        

"""
Sets the noise level of the solver CPLX_RESULT_STREAM can be::

 - *None* silent i.e. no output
 - *file* set solver to silent and output logs to *CPLX_RESULT_STREAM_FILE* cplex_output.log
 - *iostream* set solver to silent and output logs to *CPLX_RESULT_STREAM_IO* cStringIO
 - *'default'* or anything else noisy with full output closes STREAM_IO and STREAM_FILE (default)
"""    
global CPLX_RESULT_STREAM, CPLX_RESULT_STREAM_FILE, CPLX_RESULT_STREAM_IO, CPLX_SILENT_MODE
CPLX_RESULT_STREAM = 'default'
CPLX_RESULT_STREAM_FILE = None
CPLX_RESULT_STREAM_IO = None
CPLX_SILENT_MODE = True

def cplx_analyzeModel(f, lpFname=None, return_lp_obj=False, with_reduced_costs=True, with_sensitivity=False, del_intermediate=False, build_n=True, quiet=False):
    """
    Optimize a model and add the result of the optimization to the model object 
    (e.g. `reaction.value`, `objectiveFunction.value`) and the stoichiometric quiet
    matrix is automatically generated. This is a common function available in all solver interfaces.
    By default returns the model object **testing**
    
     - *f* an instantiated PySCeSCBM model object
     - *lpFname* [default=None] the name of the intermediate LP file. If not specified '_cplxtmp_.tmp'
     - *return_lp_obj* [default=False] off by default when enabled it returns the CPLEX LP object
     - *with_reduced_costs* [default=True] calculate and add reduced cost information to mode
     - *with_sensitivity* [default=False] add solution sensitivity information (not yet implemented)
     - *del_intermediate* [default=False] delete the intermediary files after updating model object, useful for server applications
     - *build_n* [default=True] generate stoichiometry from the reaction network (reactions/reagents/species)
     - *quiet* [default=False] suppress cplex output
    """

    if build_n:
        CBTools.addStoichToFBAModel(f)
    fid = f.id
    if lpFname == None:
        ##  lpFname = '_cplxtmp_%s_.tmp' % random.randint(0,100)
        lpFname = '_cplxtmp_.tmp'
    f.id = lpFname
    LPF = CBWrite.WriteModelLP(f, work_dir=os.getcwd(), quiet=quiet)
    LPF = os.path.abspath(LPF)
    flp = cplx_getModelFromLP(LPF)
    f.id = fid
    global CPLX_SILENT_MODE
    if not quiet and CPLX_SILENT_MODE:
        cplx_setOutputStreams(flp, mode='default')
    elif quiet and CPLX_SILENT_MODE:
        pass
    elif quiet:
        cplx_setOutputStreams(flp, mode=None)
    cplx_Solve(flp)
    cplx_setFBAsolutionToModel(f, flp, with_reduced_costs=with_reduced_costs)
    f.SOLUTION_STATUS_INT =  flp.solution.get_status()
    f.SOLUTION_STATUS =  cplx_getSolutionStatus(flp)
    if del_intermediate:
        os.remove(LPF)
    if return_lp_obj:
        return flp
    else:
        del flp
        return f
        
    
def cplx_setFBAsolutionToModel(fba, lp, with_reduced_costs=False):
    """
    Sets the FBA solution from a CPLEX solution to an FBA object

     - *fba* and fba object
     - *lp* a CPLEX LP object
     - *with_reduced_costs* a boolean [default=False]

    """
    sol, objname, objval = cplx_getOptimalSolution(lp)
    if not CPLX_SILENT_MODE:
        print 'Status:', cplx_getSolutionStatus(lp)
    fba.objectives[fba.activeObjIdx].solution, fba.objectives[fba.activeObjIdx].value = sol, objval
    for r in fba.reactions:
        rid = r.getPid()
        if sol.has_key(rid):
            r.value = sol[rid]
        else:
            r.value = None
    if objval != None and sol != {} and with_reduced_costs:
        RC = cplx_getReducedCosts(lp)
        CBTools.setReducedCosts(fba, RC)
    else:
        CBTools.setReducedCosts(fba, {})

def cplx_setOutputStreams(lp, mode='default'):
    """
    Sets the noise level of the solver, mode can be one of:

     - *None* silent i.e. no output
     - *'file'* set solver to silent and output logs to *CPLX_RESULT_STREAM_FILE* cplex_output.log
     - *'iostream'* set solver to silent and output logs to *CPLX_RESULT_STREAM_IO* cStringIO
     - *'default'* or anything else noisy with full output closes STREAM_IO and STREAM_FILE (default)
    """    
    global CPLX_RESULT_STREAM
    global CPLX_SILENT_MODE
    global CPLX_RESULT_STREAM_FILE
    global CPLX_RESULT_STREAM_IO
    CPLX_RESULT_STREAM = mode

    if CPLX_RESULT_STREAM == None:
        lp.set_log_stream(None)
        lp.set_results_stream(None)
        lp.set_warning_stream(None)
        CPLX_SILENT_MODE = True
    elif CPLX_RESULT_STREAM == 'file':
        if CPLX_RESULT_STREAM_FILE == None:
            CPLX_RESULT_STREAM_FILE = file('cplex_output.log', 'w')
        lp.set_log_stream(CPLX_RESULT_STREAM_FILE)
        lp.set_results_stream(CPLX_RESULT_STREAM_FILE)
        CPLX_SILENT_MODE = True
    elif CPLX_RESULT_STREAM == 'iostream':
        if CPLX_RESULT_STREAM_IO == None:
            CPLX_RESULT_STREAM_IO = cStringIO.StringIO()
        lp.set_log_stream(CPLX_RESULT_STREAM_IO)
        lp.set_results_stream(CPLX_RESULT_STREAM_IO)
        CPLX_SILENT_MODE = True
    else:
        CPLX_RESULT_STREAM = 'default'
        CPLX_SILENT_MODE = False
        try:
            CPLX_RESULT_STREAM_IO.close()
        except AttributeError:
            pass
        try:
            CPLX_RESULT_STREAM_FILE.close()
        except AttributeError:
            pass

        CPLX_RESULT_SsensitivitiesTREAM_IO = None
        CPLX_RESULT_STREAM_FILE = None
        lp.set_log_stream(os.sys.stdout)
        lp.set_results_stream(os.sys.stdout)

def cplx_getModelFromLP(lptFile, Dir=None):
    """
    Load a LPT (CPLEX format) file and return a CPLX LP model
    
     - *lptfile* an CPLEX LP format file
     - *Dir* an optional directory
    
    """
    if Dir != None:
        assert os.path.exists(Dir), '\nIncorrect path'
        lptFile = os.path.join(Dir, lptFile)
    lp = cplex.Cplex(lptFile)
    return lp

def cplx_getCPLEXModelFromLP(lptFile, Dir=None):
    """
    Load a LPT (CPLEX format) file and return a CPLX LP model
    
     - *lptfile* an CPLEX LP format file
     - *Dir* an optional directory
    
    """
    return cplx_getModelFromLP(lptFile, Dir)


def cplx_getModelFromObj(fba):
    """
    Return a GLPK object from a FBA model object (via LP file)
    """
    modId = fba.id
    ##  fba.id = '_cplxtmp_%s_.tmp' % random.randint(0,5)
    fba.id = '_cplxtmp_.tmp'
    LPF = CBWrite.WriteModelLP(fba, work_dir=os.getcwd())
    fba.id = modId
    return cplx_getModelFromLP(LPF)

def cplx_Solve(c, method='o'):
    alg = c.parameters.lpmethod.values

    if method == "o":
        c.parameters.lpmethod.set(alg.auto)
    elif method == "p":
        c.parameters.lpmethod.set(alg.primal)
    elif method == "d":
        c.parameters.lpmethod.set(alg.dual)
    elif method == "b":
        c.parameters.lpmethod.set(alg.barrier)
        c.parameters.barrier.crossover.set(c.parameters.barrier.crossover.values.none)
    elif method == "h":
        c.parameters.lpmethod.set(alg.barrier)
    elif method == "s":
        c.parameters.lpmethod.set(alg.sifting)
    elif method == "c":
        c.parameters.lpmethod.set(alg.concurrent)
    else:
        print "Unrecognized option, using automatic"
        c.parameters.lpmethod.set(alg.auto)
    try:
        c.solve()
    except cplex.exceptions.CplexSolverError:
        print "Exception raised during solve"
        return None
    # solution.get_status() returns an integer code
    status = c.solution.get_status()
    if status == c.solution.status.optimal:
        if not CPLX_SILENT_MODE:
            print "INFO: Model is optimal:", status
    elif status == c.solution.status.unbounded:
        if not CPLX_SILENT_MODE:
            print "INFO: Model is unbounded:", status
        return status
    elif status == c.solution.status.infeasible:
        if not CPLX_SILENT_MODE:        
            print "INFO: Model is infeasible:", status
        return status
    elif status == c.solution.status.infeasible_or_unbounded:
        if not CPLX_SILENT_MODE:        
            print "INFO: Model is infeasible or unbounded:", status
        return status    

    if not CPLX_SILENT_MODE:
        s_method = c.solution.get_method()
        print "Solution status = " , status, ":",
        # the following line prints the status as a string
        print c.solution.status[status]
        print "Solution method = ", s_method, ":",
        print c.solution.method[s_method]
    
    s_type = c.solution.get_solution_type()
    if s_type == c.solution.type.none:
        print "INFO: No solution available"
        return None
    else:
        if not CPLX_SILENT_MODE:
            print "Objective value = " , c.solution.get_objective_value()
        return status

def cplx_SolveMILP(c, method=None):
    try:
        c.solve()
    except cplex.exceptions.CplexSolverError:
        print "Exception raised during solve"

def cplx_getOptimalSolution(c):
    """
    From a CPLX model extract a tuple of solution, ObjFuncName and ObjFuncVal
    """
    s_val = []
    s_name = []
    fba_sol = {}
    objf_name = None
    objf_val = None
    try:
        objf_name = c.objective.get_name()        
        objf_val = c.solution.get_objective_value()        
        s_val = c.solution.get_values()
        s_name = c.variables.get_names()
        for n in range(len(s_name)):
            fba_sol.update({s_name[n] : s_val[n]})        
    except Exception, ex:
        print ex
        print 'WARNING: No solution to get'
        s_val = []
        s_name = []
        fba_sol = {}
        objf_val = None
    del s_name,s_val
    return fba_sol, objf_name, objf_val

def cplx_getOptimalSolution2(c, names):
    """
    From a CPLX model extract a tuple of solution, ObjFuncName and ObjFuncVal
    """
    s_val = []
    s_name = []
    fba_sol = None
    objf_name = None
    objf_val = None
    try:
        s_val = c.solution.get_values()
        objf_name = c.objective.get_name()
        objf_val = c.solution.get_objective_value()        
        s_name = names
    except AttributeError:
        print 'No solution attribute on CPLEX object'    
    fba_sol = {}
    for n in range(len(s_name)):
        fba_sol.update({s_name[n] : s_val[n]})
    del names,s_val
    return fba_sol, objf_name, objf_val

def cplx_writeLPsolution(fba_sol, objf_name, objf_val, fname, Dir=None, separator=','):
    """
    This function writes the optimal solution, produced wth `cplx_getOptimalSolution` to file
    
     - *fba_sol* a dictionary of Flux : value pairs
     - *objf_name* the objective function flux id
     - *objf_val* the optimal value of the objective function
     - *fname* the output filename
     - *Dir* [default=None] use directory if not None
     - *separator* [default=','] the column separator
     
    """
    if Dir != None:
        assert os.path.exists(Dir), '\nPath does not exist'
        fname = os.path.join(Dir, fname)
    fname += '.csv'
    F = file(fname, 'w')
    cntr = 0
    F.write('%s%s%s\n' % ('ObjectiveFunction', separator, objf_name))
    for r in fba_sol:
        F.write('%s%s%f\n' % (r, separator, fba_sol[r]))
    F.flush()
    F.close()
    print 'CSV exported to %s' % fname

def cplx_getReducedCosts(c):
    """
    Extract ReducedCosts from LP and return as a dictionary

    - CPLEX LP

    """
    s_name = c.variables.get_names()
    r_costs = c.solution.get_reduced_costs()
    output = {}
    for v in range(len(s_name)):
        output.update({s_name[v] : r_costs[v]})
    del s_name, r_costs
    return output

def cplx_getSensitivities(c):
    """
    Get the sensitivities of each constraint on the objective function with input

     - *c* a CPLEX LP

    Output is a tuple of bound and objective sensitivities where the objective
    sensitivity is described in the CPLEX reference manual as::

     ... the objective sensitivity shows each variable, its reduced cost and the range over
     which its objective function coefficient can vary without forcing a change
     in the optimal basis. The current value of each objective coefficient is
     also displayed for reference.

     - *objective coefficient sensitivity* {flux : (reduced_cost, lower_obj_sensitivity, coeff_value, upper_obj_sensitivity)}
     - *rhs sensitivity* {constraint : (low, value, high)}
     - *bound sensitivity ranges* {flux : (lb_low, lb_high, ub_low, ub_high)}

    """
    SENSE_RHS = {}
    SENSE_BND = {}
    SENSE_OBJ = {}
    c_names = c.linear_constraints.get_names()
    rhs_val = c.linear_constraints.get_rhs()
    j_names = c.variables.get_names()

    rhs_sense = c.solution.sensitivity.rhs()
    bnd_sense = c.solution.sensitivity.bounds()
    obj_sense = c.solution.sensitivity.objective()
    obj_coeff = c.objective.get_linear()
    red_cost = c.solution.get_reduced_costs()

    for r in range(c.variables.get_num()):
        SENSE_BND.update({j_names[r] : (bnd_sense[r][0], bnd_sense[r][1], bnd_sense[r][2], bnd_sense[r][3])})
        SENSE_OBJ.update({j_names[r] : (red_cost[r], obj_sense[r][0], obj_coeff[r], obj_sense[r][1])})

    for s in range(c.linear_constraints.get_num()):
        SENSE_RHS.update({c_names[s] : (rhs_sense[s][0], rhs_val[s], rhs_sense[s][1])})

    return (SENSE_OBJ, SENSE_RHS, SENSE_BND)

def cplx_getDualValues(c):
    """
    Get the get the dual values of the solution

     - *c* a CPLEX LP

    Output is a dictionary of {name : value} pairs

    """
    d_names = c.linear_constraints.get_names()
    d_values = c.solution.get_dual_values()
    output = {}
    for j in range(len(d_names)):
        output.update({d_names[j] : d_values[j]})
    return output


def cplx_getSolutionStatus(c):
    """
    Returns one of:

     - *LPS_OPT*: solution is optimal;
     - *LPS_FEAS*: solution is feasible;
     - *LPS_INFEAS*: solution is infeasible;
     - *LPS_NOFEAS*: problem has no feasible solution;
     - *LPS_UNBND*: problem has unbounded solution;
     - *LPS_UNDEF*: solution is undefined.
     - *LPS_NONE*: no solution

    """
    # solution.get_status() returns an integer code
    if c.solution.get_solution_type() == c.solution.type.none:
        print "\nNo solution available\n"
        return 'LPS_NONE'
    status = c.solution.get_status()
    if status == c.solution.status.optimal:
        if not CPLX_SILENT_MODE:
            print "Model is optimal"
        return 'LPS_OPT'
    elif status == c.solution.status.feasible:
        print "\nModel is feasible"
        return 'LPS_FEAS'
    elif status == c.solution.status.unbounded:
        print "\nModel is unbounded"
        return 'LPS_UNBND'
    elif status == c.solution.status.infeasible:
        print "\nModel is infeasible"
        return 'LPS_INFEAS'
    elif status == c.solution.status.infeasible_or_unbounded:
        print "\nModel is infeasible or unbounded"
        return 'LPS_INFEAS or LPS_UNBND'
    else:
        return 'LPS_NONE'

def cplx_setObjective(c, id, expr=None, sense='maximize', reset=True):
    """
    Set a new objective function note that there is a major memory leak in
    `c.variables.get_names()` whch is used when reset=True. If this is a problem
    use cplx_setObjective2 which takes *names* as an input:
    
     - *c* a CPLEX LP object
     - *id* the r_id of the flux to be optimized
     - *expr* a list of (coefficient, flux) pairs
     - *sense* 'maximize'/'minimize'
     - *reset* [default=True] reset all objective function coefficients to zero
    
    """
    sense = sense.lower()
    if sense == 'max': sense = 'maximize'
    if sense == 'min': sense = 'minimize'
    if sense in ['maximise', 'minimise']:
        sense = sense.replace('se','ze')
    assert sense in ['maximize', 'minimize'], "\nsense must be ['maximize', 'minimize'] not %s" % sense    
    c.objective.set_name(id)
    if reset:
        vars = c.variables.get_names()
        new_obj = []
        for e in range(c.variables.get_num()):
            new_obj.append((vars[e], 0.0))
        c.objective.set_linear(new_obj)
        del new_obj, vars
    if expr != None:
        expr = [(e[1], e[0]) for e in expr]
        if __DEBUG__: print expr
        c.objective.set_linear(expr)
        del expr

    if sense == 'minimize':
        c.objective.set_sense(c.objective.sense.minimize)
        if __DEBUG__: print 'Set minimizing'
    else:
        c.objective.set_sense(c.objective.sense.maximize)
        if __DEBUG__: print 'Set maximizing'

def cplx_setObjective2(c, id, names, expr=None, sense='maximize', reset=True):
    """
    Set a new objective function. This is a workaround function to avoid the
    e is a major memory leak in `c.variables.get_names()` whch is used
    in cplx_setObjective()  when reset=True.
    
    """
    sense = sense.lower()
    if sense == 'max': sense = 'maximize'
    if sense == 'min': sense = 'minimize'    
    if sense in ['maximise', 'minimise']:
        sense = sense.replace('se','ze')
    assert sense in ['maximize', 'minimize'], "\nsense must be ['maximize', 'minimize'] not %s" % sense        
    c.objective.set_name(id)
    if reset:
        ##  var_num = c.variables.get_num()
        ##  names = range(var_num)
        new_obj = []
        for e in range(len(names)):
            new_obj.append((names[e], 0.0))
        c.objective.set_linear(new_obj)
        del new_obj, names
    if expr != None:
        expr = [(e[1], e[0]) for e in expr]
        if __DEBUG__: print expr
        c.objective.set_linear(expr)
        del expr
    if sense == 'minimize':
        c.objective.set_sense(c.objective.sense.minimize)
        if __DEBUG__: print 'Set minimizing'
    else:
        c.objective.set_sense(c.objective.sense.maximize)
        if __DEBUG__: print 'Set maximizing'


def cplx_setSingleConstraint(c, id, expr=[], sense='E', rhs=0.0):
    ind = []
    val = []
    for e in expr:
        ind.append(e[1])
        val.append(e[0])
    lexp = cplex.SparsePair(ind, val)
    c.linear_constraints.add(lin_expr=[lexp], senses=sense, rhs=[rhs], names=[id])
    del ind, val, lexp

def cplx_setBounds(c, id, min=None, max=None):
    if min != None:
        c.variables.set_lower_bounds(id, min)
    if max != None:
        c.variables.set_upper_bounds(id, max)

def cplx_writeLPtoLPTfile(c, filename, title=None, Dir=None):
    """
    Write out a CPLEX model as an LP format file
    
    """
    if Dir != None:
        filename = os.path.join(Dir, filename)
    if title != None:
        c.set_problem_name(title)
    c.write(filename+'.lp', filetype='lp')
    print 'LP output as %s' % filename+'.lp'

def cplx_runInputScan(fba, exDict, wDir, input_lb=-10.0, input_ub=0.0, writeHformat=False, rationalLPout=False):
    """
    scans all inputs
    
    """
    debug_dir = os.path.join(wDir, 'lpt')
    ine_dir = os.path.join(wDir, 'ine')
    rat_dir = os.path.join(wDir, 'eslv')
    if not os.path.exists(debug_dir):
        os.mkdir(debug_dir)
    if writeHformat and not os.path.exists(ine_dir):
        os.mkdir(ine_dir)
    if rationalLPout and not os.path.exists(rat_dir):
        os.mkdir(rat_dir)

    optimal_growth_rates = {}
    infeasible_inputs = []
    modname = fba.inputfile_id
    fbaid0 = fba.id
    CBTools.addStoichToFBAModel(fba)
    for input in exDict.keys():
        fba.id = '%s(%s)' % (modname, input)
        ilb = fba.getFluxBoundByReactionID(input, 'lower').value
        iub = fba.getFluxBoundByReactionID(input, 'upper').value
        fba.setReactionBound(input, input_lb, 'lower')
        fba.setReactionBound(input, input_ub, 'upper')
        fbalp = cplx_getModelFromObj(fba)
        cplx_Solve(fbalp)
        if cplx_getSolutionStatus(fbalp) == 'LPS_OPT':
            opt = cplx_getOptimalSolution(fbalp)[2]
            if __DEBUG__: print 'Optimal growth rate %s(%s) = %s' % (modname, input, opt)
            optimal_growth_rates.update({input : opt})
            if writeHformat:
                CBWrite.WriteModelHFormatFBA(fba, ine_dir, realnum=False)
                CBWrite.WriteModelHFormatFBA(fba, ine_dir, realnum=True)
            if rationalLPout:
                tmp = CBWrite.WriteModelLP(fba, work_dir=rat_dir)
        else:
            print 'CPLX returned an error code: %s' % cplx_getSolutionStatus(fbalp)
            infeasible_inputs.append(input)
        fba.setReactionBound(input, ilb, 'lower')
        fba.setReactionBound(input, iub, 'upper')
        del fbalp
    fba.id = fbaid0
    return optimal_growth_rates, infeasible_inputs

def cplx_runInputScanWithFVA(fba, exDict, wDir, input_lb=-10.0, input_ub=0.0, writeHformat=False, rationalLPout=False, NaNonFail=False):
    """
    scans all inputs
    
    """
    raise NotImplementedError, 'Not yet anyway'

    debug_dir = os.path.join(wDir, 'lpt')
    ine_dir = os.path.join(wDir, 'ine')
    rat_dir = os.path.join(wDir, 'eslv')
    if not os.path.exists(debug_dir):
        os.mkdir(debug_dir)
    if writeHformat and not os.path.exists(ine_dir):
        os.mkdir(ine_dir)
    if rationalLPout and not os.path.exists(rat_dir):
        os.mkdir(rat_dir)

    optimal_growth_rates = {}
    infeasible_inputs = []
    modname = fba.inputfile_id
    fbaid0 = fba.id
    fva_out = []
    for input in exDict.keys():
        fba.id = '%s(%s)' % (modname, input)
        ilb = fba.getFluxBoundByReactionID(input, 'lower').value
        iub = fba.getFluxBoundByReactionID(input, 'upper').value
        fba.setReactionBound(input, input_lb, 'lower')
        fba.setReactionBound(input, input_ub, 'upper')
        LPfile = CBWrite.WriteModelLP(fba, work_dir=debug_dir)
        fbalp = PyscesGLPK.glpk_getModelFromLPT(LPfile, debug_dir)
        ##  fbalp.solve()
        PyscesGLPK.glpk_Solve(fbalp)
        if fbalp.solution_status == 'LPS_OPT':
            opt = PyscesGLPK.glpk_getOptimalSolution(fbalp)[2]
            if __DEBUG__: print 'Optimal growth rate %s(%s) = %s' % (modname, input, opt)
            optimal_growth_rates.update({input : opt})
            if writeHformat:
                CBWrite.WriteModelHFormatFBA(fba, ine_dir, realnum=False)
                CBWrite.WriteModelHFormatFBA(fba, ine_dir, realnum=True)
            if rationalLPout:
                tmp = CBWrite.WriteModelLP(fba, work_dir=rat_dir)
            PyscesGLPK.glpk_addOptimalSolutionAsConstraint(fbalp, tol=None, bound='lower', optPercentage=100.0, reset_of=True)
            OUTPUT_ARRAY, ncol_names = doFluxVariability(fbalp, ignore_fluxes=[], write_debug_lps=False,\
                                        pause_on_skip=True, model_name=fba.id, work_dir=wDir)
            fva_out.append((input,OUTPUT_ARRAY, ncol_names))
            del OUTPUT_ARRAY, ncol_names
        else:
            print 'GLPK returned an error code: %s' % fbalp.solution_status
            infeasible_inputs.append(input)
            os.remove(LPfile)
        fba.setReactionBound(input, ilb, 'lower')
        fba.setReactionBound(input, iub, 'upper')
    fba.id = fbaid0
    return optimal_growth_rates, infeasible_inputs, fva_out

def cplx_func_GetCPXandPresolve(fba, pre_opt, objF2constr, quiet=False):
    """
    This is a utility function that does a presolve for FVA, MSAF etc. Generates properly formatted
    empty objects if pre_opt == False
    
     - *pre_opt* a boolean 
     - *fba* a CBModel object
     - *objF2constr* add objective function as constraint
     - *quiet* [default=False] supress cplex output
     
    Returns: pre_sol, pre_oid, pre_oval, OPTIMAL_PRESOLUTION, REDUCED_COSTS
     
    """
    cpx = cplx_getModelFromObj(fba)
    if quiet:
        cplx_setOutputStreams(cpx, mode=None)    
    OPTIMAL_PRESOLUTION = None
    pre_sol = pre_oid = pre_oval = None
    REDUCED_COSTS = {}
    if pre_opt:
        status = cplx_Solve(cpx, method='o')
        if cplx_getSolutionStatus(cpx) == 'LPS_OPT':
            print 'Valid Presolution'
            OPTIMAL_PRESOLUTION = True
            pre_sol, pre_oid, pre_oval = cplx_getOptimalSolution(cpx)
            fba.objectives[fba.activeObjIdx].solution, fba.objectives[fba.activeObjIdx].value = pre_sol, pre_oval
            REDUCED_COSTS = cplx_getReducedCosts(cpx)
        else:
            print 'Invalid Presolution'
            OPTIMAL_PRESOLUTION = False
            pre_sol = {}
            for r in fba.reactions:
                pre_sol.update({r.getPid : 0.0})
                r.reduced_cost = 0.0
            pre_oval = 0.0
            pre_oid = 'None'
            raise RuntimeError, '\nPresolve failed to optimize this model and cannot continue!'
    else:
        pre_sol = {}
        for r in fba.reactions:
            pre_sol.update({r.getPid() : 0.0})
        if objF2constr:
            pre_oval = fba.objectives[fba.activeObjIdx].value
            pre_oid = fba.objectives[fba.activeObjIdx].getPid()
        else:
            pre_oval = 0.0
            pre_oid = 'None'
    for r in fba.reactions:
        r.value = pre_sol[r.id]
    return cpx, pre_sol, pre_oid, pre_oval, OPTIMAL_PRESOLUTION, REDUCED_COSTS

def cplx_func_SetObjectiveFunctionAsConstraint(cpx, rhs_sense, oval, tol, optPercentage):
    """
    Take the objective function and "optimum" value and add it as a constraint
     - *cpx* a cplex object
     - *oval* the objective value
     - *tol*  [default=None] do not floor/ceiling the objective function constraint, otherwise round of to *tol*
     - *rhs_sense* [default='lower'] means objC >= objVal the inequality to use for the objective constraint can also be *upper* or *equal*
     - *optPercentage* [default=100.0] means the percentage optimal value to use for the RHS of the objective constraint: optimal_value*(optPercentage/100.0)
    
    """
    
    # generate new constraint from old objective value (use non-zero coefficients)
    LCS = cpx.objective.get_linear()
    LCN = cpx.variables.get_names()
    new_constraint = []
    for c in range(len(LCS)):
        if LCS[c] != 0.0:
            new_constraint.append((LCS[c], LCN[c]))

    # set objective constraint
    if rhs_sense == 'equal':
        cplx_setSingleConstraint(cpx, 'ObjCstr', expr=new_constraint, sense='E', rhs=oval)
        ##  cplx_setSingleConstraint(cpx, 'ObjCstr', expr=[(1, pre_oid)], sense='E', rhs=oval)
    elif rhs_sense == 'upper':
        if tol != None:
            ub =  numpy.ceil(oval/tol)*tol*optPercentage/100.0
        else:
            ub = oval*(optPercentage/100.0)
        cplx_setSingleConstraint(cpx, 'ObjCstr', expr=new_constraint, sense='L', rhs=ub)
        ##  cplx_setSingleConstraint(cpx, 'ObjCstr', expr=[(1, pre_oid)], sense='L', rhs=ub)
    elif rhs_sense == 'lower':
        if tol != None:
            lb =  numpy.floor(oval/tol)*tol*optPercentage/100.0
        else:
            lb = oval*(optPercentage/100.0)
        cplx_setSingleConstraint(cpx, 'ObjCstr', expr=new_constraint, sense='G', rhs=lb)
        ##  cplx_setSingleConstraint(cpx, 'ObjCstr', expr=[(1, pre_oid)], sense='G', rhs=lb)
    else:
        raise RuntimeError, "\nInvalid RHS sense: %s" % rhs_sense

def cplx_MinimizeSumOfAbsFluxes(fba, selected_reactions=None, pre_opt=True, tol=None, objF2constr=True, rhs_sense='lower', optPercentage=100.0, work_dir=None, quiet=False, debug=False, objective_coefficients={}, return_lp_obj=False):
    """
    Minimize the sum of absolute fluxes sum(abs(J1) + abs(J2) + abs(J3) ... abs(Jn)) by adding two constraints per flux
    and a variable representing the absolute value:
 
     Min: Ci abs_Ji
      Ji - abs_Ji <= 0
      Ji + abs_Ji >= 0 
      
     Such that:
      NJi = 0
      Jopt = opt

    Arguments:

     - *fba* an FBA model object
     - *selected reactions* [default=None] means use all reactions otherwise use the reactions listed here
     - *pre_opt* [default=True] attempt to presolve the FBA and report its results in the ouput, if this is diabled and *objF2constr* is True then the vid/value of the current active objective is used
     - *tol*  [default=None] do not floor/ceiling the objective function constraint, otherwise round of to *tol*
     - *rhs_sense* [default='lower'] means objC >= objVal the inequality to use for the objective constraint can also be *upper* or *equal*
     - *optPercentage* [default=100.0] means the percentage optimal value to use for the RHS of the objective constraint: optimal_value*(optPercentage/100.0)
     - *work_dir* [default=None] the MSAF working directory for temporary files default = cwd+fva
     - *debug* [default=False] if True write out all the intermediate MSAF LP's into work_dir
     - *quiet* [default=False] if enabled supress CPLEX output
     - *objF2constr* [default=True] add the model objective function as a constraint using rhs_sense etc. If
       this is True with pre_opt=False then the id/value of the active objective is used to form the constraint
     - *objective_coefficients* [default={}] a dictionary of (reaction_id : float) pairs that provide the are introduced as objective coefficients to the absolute flux value. Note that the default value of the coefficient (non-specified) is +1.
     - *return_lp_obj* [default=False] off by default when enabled it returns the CPLEX LP object    

    With outputs:
    
     - *fba* an update instance of a CBModel. Note that the FBA model objective function value is the original value set as a constraint
    
    """
    if work_dir == None:
        work_dir = os.getcwd()
    else:
        assert os.path.exists(work_dir), '\nWhat did you think would happen now!'
    if debug:
        debug_dir = os.path.join(work_dir,'DEBUG')
        if not os.path.exists(debug_dir):
            os.mkdir(debug_dir)
            
    # generate a presolution
    cpx = OPTIMAL_PRESOLUTION = None
    pre_sol = pre_oid = pre_oval = None
    REDUCED_COSTS = {}
    cpx, pre_sol, pre_oid, pre_oval, OPTIMAL_PRESOLUTION, REDUCED_COSTS = cplx_func_GetCPXandPresolve(fba, pre_opt, objF2constr, quiet=quiet)
    # if required add the objective function as a constraint
    if objF2constr:
        cplx_func_SetObjectiveFunctionAsConstraint(cpx, rhs_sense, pre_oval, tol, optPercentage)

    STORED_OPT = None
    if pre_opt:
        STORED_OPT = pre_oval
    else:
        STORED_OPT = fba.getActiveObjective().getValue()

    # minimize the absolute sum of fluxes
    if selected_reactions != None:
        rids = fba.getReactionIds()
        for r in selected_reactions:
            assert r in rids, "\n%s is not a valid reaction name" % r
    else:
        selected_reactions = fba.getReactionIds()

    # this removes the objective function 
    ##  fba_obj_ids = fba.getActiveObjective().getFluxObjectiveReactions()
    ##  print fba_obj_ids
    ##  for R in fba_obj_ids:
        ##  if R in selected_reactions:
            ##  selected_reactions.pop(selected_reactions.index(R))
    ##  del R, fba_obj_ids
    NUM_FLX = len(selected_reactions)
    print 'Total number of reactions: %s' % NUM_FLX
    abs_selected_reactions = ['abs_%s' % r for r in selected_reactions]
    # absJ implicitly assumed to be absJ >= 0
    ##  lbs = [0.0]*len(abs_selected_reactions)
    ##  ubs = [numpy.inf]*len(abs_selected_reactions)
    ##  cpx.variables.add(lb=lbs, ub=ubs, names=abs_selected_reactions)
    cpx.variables.add(names=abs_selected_reactions)
    """
    cpx.variables.add(obj=[], lb=[], ub=[], types='', names=[], columns=[])
    """
    ##  obj_func = [(1, r) for r in abs_selected_reactions]

    """
     J - abs_J <= 0
     J + abs_J >= 0     
    """
    lin_exprs = []
    senses = []
    names = []
    rhss = []
    range_values = []  
    obj_func = []
    
    for r in range(len(selected_reactions)):
        name1 = abs_selected_reactions[r]+'_c1'
        lin_expr1 = cplex.SparsePair(ind = [selected_reactions[r], abs_selected_reactions[r]], val = [1.0, -1.0])
        sense1 = 'L'
        name2 = abs_selected_reactions[r]+'_c2'
        lin_expr2 = cplex.SparsePair(ind = [selected_reactions[r], abs_selected_reactions[r]], val = [1.0, 1.0])
        sense2 = 'G'
        names.append(name1)
        lin_exprs.append(lin_expr1)
        senses.append(sense1)
        rhss.append(0.0)
        range_values.append(0.0)
        names.append(name2)
        lin_exprs.append(lin_expr2)
        senses.append(sense2)
        rhss.append(0.0)
        range_values.append(0.0)
        ##  print selected_reactions[r], abs_selected_reactions[r]
        if objective_coefficients.has_key(selected_reactions[r]):
            obj_func.append((objective_coefficients[selected_reactions[r]], abs_selected_reactions[r]))
        else:
            obj_func.append((1, abs_selected_reactions[r]))
            
    cplx_setObjective(cpx, 'MAFS', obj_func, 'min', reset=True)  
    cpx.linear_constraints.add(lin_expr=lin_exprs, senses=senses, rhs=rhss, range_values=range_values, names=names)    
    if debug:
        cplx_writeLPtoLPTfile(cpx, 'MSAF_base_(%s)' % time.time(), title=None, Dir=debug_dir)
        ##  cplx_writeLPtoLPTfile(cpx, 'MSAF_base_%s' % time.time() , title=None, Dir=debug_dir)
    
    cplx_Solve(cpx, method='o')
    cplx_setFBAsolutionToModel(fba, cpx, with_reduced_costs=True)
    
    print '\nFlux minimization objective value:', fba.getActiveObjective().getValue()
    fba.getActiveObjective().setValue(STORED_OPT)
    if quiet:
        cplx_setOutputStreams(cpx, mode='default')
    if not return_lp_obj:
        return fba
    else:
        return cpx


def cplx_MinimizeNumActiveFluxes(fba, selected_reactions=None, pre_opt=True, tol=None, objF2constr=True, rhs_sense='lower', optPercentage=100.0, work_dir=None, quiet=False, debug=False, objective_coefficients={}, return_lp_obj=False, populate=None):
    """
    Minimize the sum of active fluxes
 
     Min: sum(Bi)
      Bi = 0 -> Ci Ji = 0
      
     Such that:
      NJi = 0
      Jbio = opt
      
     where:
       Binary Bi    

    Arguments:

     - *fba* an FBA model object
     - *selected reactions* [default=None] means use all reactions otherwise use the reactions listed here
     - *pre_opt* [default=True] attempt to presolve the FBA and report its results in the ouput, if this is diabled and *objF2constr* is True then the vid/value of the current active objective is used
     - *tol*  [default=None] do not floor/ceiling the objective function constraint, otherwise round of to *tol*
     - *rhs_sense* [default='lower'] means objC >= objVal the inequality to use for the objective constraint can also be *upper* or *equal*
     - *optPercentage* [default=100.0] means the percentage optimal value to use for the RHS of the objective constraint: optimal_value*(optPercentage/100.0)
     - *work_dir* [default=None] the MSAF working directory for temporary files default = cwd+fva
     - *debug* [default=False] if True write out all the intermediate MSAF LP's into work_dir
     - *quiet* [default=False] if enabled supress CPLEX output
     - *objF2constr* [default=True] add the model objective function as a constraint using rhs_sense etc. If
       this is True with pre_opt=False then the id/value of the active objective is used to form the constraint
     - *objective_coefficients* [default={}] a dictionary of (reaction_id : float) pairs that provide the are introduced as objective coefficients to the absolute flux value. Note that the default value of the coefficient (non-specified) is +1.
     - *return_lp_obj* [default=False] off by default when enabled it returns the CPLEX LP object
     - *populate* [default=None] enable search algorithm to find multiple (sub)optimal solutions. Set with a tuple of (RELGAP=0.0, POPULATE_LIMIT=20, TIME_LIMIT=300) suggested values only.
      - *RELGAP* [default=0.0] relative gap to optimal solution
      - *POPULATE_LIMIT * [default=20] terminate when so many solutions have been found
      - *TIME_LIMIT* [default=300] terminate search after so many seconds important with higher values of "POPULATION_LIMIT*
    
    With outputs:
    
     - *fba* an update instance of a CBModel OR
     - *cpx* a cplex model OR
     - *populate_data* a population data set OR
     - *cpx, populate_data* both the cps object and population data set
     
    depending on selected flags.
    
    """   
    
    if work_dir == None:
        work_dir = os.getcwd()
    else:
        assert os.path.exists(work_dir), '\nWhat did you think would happen now!'
    if debug:
        debug_dir = os.path.join(work_dir,'DEBUG')
        if not os.path.exists(debug_dir):
            os.mkdir(debug_dir)
            
    # generate a presolution
    cpx = OPTIMAL_PRESOLUTION = None
    pre_sol = pre_oid = pre_oval = None
    REDUCED_COSTS = {}
    cpx, pre_sol, pre_oid, pre_oval, OPTIMAL_PRESOLUTION, REDUCED_COSTS = cplx_func_GetCPXandPresolve(fba, pre_opt, objF2constr, quiet=quiet)
    # if required add the objective function as a constraint
    if objF2constr:
        cplx_func_SetObjectiveFunctionAsConstraint(cpx, rhs_sense, pre_oval, tol, optPercentage)

    # minimize the absolute sum of fluxes
    if selected_reactions != None:
        rids = fba.getReactionIds()
        for r in selected_reactions:
            assert r in rids, "\n%s is not a valid reaction name" % r
    else:
        selected_reactions = fba.getReactionIds()

    # this removes the objective function 
    ##  fba_obj_ids = fba.getActiveObjective().getFluxObjectiveReactions()
    ##  print fba_obj_ids
    ##  for R in fba_obj_ids:
        ##  if R in selected_reactions:
            ##  selected_reactions.pop(selected_reactions.index(R))
    ##  del R, fba_obj_ids
    NUM_FLX = len(selected_reactions)
    print 'Total number of reactions: %s' % NUM_FLX
    bin_selected_reactions = ['bin_%s' % r for r in selected_reactions]
    # bin_J in [0,1]
    ##  lbs = [0.0]*len(bin_selected_reactions)
    ##  ubs = [numpy.inf]*len(bin_selected_reactions)
    ##  cpx.variables.add(lb=lbs, ub=ubs, names=bin_selected_reactions)
    cpx.variables.add(names=bin_selected_reactions)
    """
    cpx.variables.add(obj=[], lb=[], ub=[], types='', names=[], columns=[])
    """
    
    for bv in bin_selected_reactions:
        cpx.variables.set_types(bv, cpx.variables.type.binary)
        
    obj_func = []
    
    for r in range(len(selected_reactions)):
        name1 = 'ind_'+selected_reactions[r]
        ##  print selected_reactions[r]
        lin_expr1 = cplex.SparsePair(ind = [selected_reactions[r]], val = [1.0])
        sense1 = 'E'
        rhs1 = 0.0
        indvar1 = cpx.variables.get_indices(bin_selected_reactions[r])
        complemented1 = 1
        cpx.indicator_constraints.add(lin_expr=lin_expr1, sense=sense1, rhs=rhs1, indvar=indvar1, complemented=complemented1, name=name1)
        obj_func.append((1, bin_selected_reactions[r]))

        if objective_coefficients.has_key(selected_reactions[r]):
            obj_func.append((objective_coefficients[selected_reactions[r]], bin_selected_reactions[r]))
        else:
            obj_func.append((1, bin_selected_reactions[r]))


    cplx_setObjective(cpx, 'MNAF', obj_func, 'min', reset=True) 
    
    if debug:
        cplx_writeLPtoLPTfile(cpx, 'MNAF_base_(%s)' % time.time() , title=None, Dir=debug_dir)
        ##  cplx_writeLPtoLPTfile(cpx, 'MNAV_base_%s' % time.time() , title=None, Dir=debug_dir)
    
    if populate == None:
        cplx_SolveMILP(cpx) #  cpx.solve()
        cplx_setFBAsolutionToModel(fba, cpx, with_reduced_costs=False)
        if quiet:
            cplx_setOutputStreams(cpx, mode='default')
        if not return_lp_obj:
            return fba
        else:
            return cpx        
    else:
        RELGAP = populate[0] #0.0 # relative gap to optimal solution
        POPULATE_LIMIT = populate[1] #20
        TIME_LIMIT =  populate[2] #300 # seconds
        INTENSITY = cpx.parameters.mip.pool.intensity.values.very_aggressive
        # ABSGAP = 0.0
        cpx.parameters.mip.pool.relgap.set(RELGAP) # Gunnar
        # cpx.parameters.mip.pool.absgap.set(ABSGAP) # check this
        cpx.parameters.mip.pool.intensity.set(INTENSITY) # get "all" (sub)optimal solutions
        cpx.parameters.mip.limits.populate.set(POPULATE_LIMIT)
        cpx.parameters.timelimit.set(TIME_LIMIT)
        cpx.populate_solution_pool()
        
        population = []
        var_names = cpx.variables.get_names()
        var_num = cpx.variables.get_num()
        pop_names = cpx.solution.pool.get_names()
        pop_num = cpx.solution.pool.get_num()
        
        population.append(var_names)
        for p in range(pop_num):
            population.append(cpx.solution.pool.get_values(p))

        if quiet:
            cplx_setOutputStreams(cpx, mode='default')
        if not return_lp_obj:
            return population
        else:
            return cpx, population 
        

def cplx_FluxVariabilityAnalysis(fba, selected_reactions=None, pre_opt=True, tol=None, objF2constr=True, rhs_sense='lower', optPercentage=100.0, work_dir=None, quiet=False, debug=False):
    """
    Perform a flux variability analysis on an fba model:

     - *fba* an FBA model object
     - *selected reactions* [default=None] means use all reactions otherwise use the reactions listed here
     - *pre_opt* [default=True] attempt to presolve the FBA and report its results in the ouput, if this is diabled and *objF2constr* is True then the vid/value of the current active objective is used
     - *tol*  [default=None] do not floor/ceiling the objective function constraint, otherwise round of to *tol*
     - *rhs_sense* [default='lower'] means objC >= objVal the inequality to use for the objective constraint can also be *upper* or *equal*
     - *optPercentage* [default=100.0] means the percentage optimal value to use for the RHS of the objective constraint: optimal_value*(optPercentage/100.0)
     - *work_dir* [default=None] the FVA working directory for temporary files default = cwd+fva
     - *debug* [default=False] if True write out all the intermediate FVA LP's into work_dir
     - *quiet* [default=False] if enabled supress CPLEX output
     - *objF2constr* [default=True] add the model objective function as a constraint using rhs_sense etc. If
       this is True with pre_opt=False then the id/value of the active objective is used to form the constraint

    and returns an array with columns::

     Reaction, Reduced Costs, Variability Min, Variability Max, abs(Max-Min), MinStatus, MaxStatus

    and a list containing the row names.
    """
    if work_dir == None:
        work_dir = os.getcwd()
    else:
        assert os.path.exists(work_dir), '\nWhat did you think would happen now!'
    if debug:
        debug_dir = os.path.join(work_dir,'DEBUG')
        if not os.path.exists(debug_dir):
            os.mkdir(debug_dir)
    
    # generate a presolution
    cpx = OPTIMAL_PRESOLUTION = None
    pre_sol = pre_oid = pre_oval = None
    REDUCED_COSTS = {}
    cpx, pre_sol, pre_oid, pre_oval, OPTIMAL_PRESOLUTION, REDUCED_COSTS = cplx_func_GetCPXandPresolve(fba, pre_opt, objF2constr, quiet=quiet)
    # if required add the objective function as a constraint
    if objF2constr:
        cplx_func_SetObjectiveFunctionAsConstraint(cpx, rhs_sense, pre_oval, tol, optPercentage)
    if debug:
        cplx_writeLPtoLPTfile(cpx, 'FVA_base', title=None, Dir=debug_dir)

    # do the FVA
    NUM_FLX = len(fba.reactions)
    print 'Total number of reactions: %s' % NUM_FLX
    if selected_reactions != None:
        rids = fba.getReactionIds()
        for r in selected_reactions:
            assert r in rids, "\n%s is not a valid reaction name" % r
    else:
        selected_reactions = fba.getReactionIds()
    NUM_FLX = len(selected_reactions)        
    print 'Number of user selected variables: %s' % NUM_FLX
    try:
        OUTPUT_ARRAY = numpy.zeros((NUM_FLX, 7), numpy.double)
    except AttributeError:
        OUTPUT_ARRAY = numpy.zeros((NUM_FLX, 7))
    OUTPUT_NAMES = []
    cntr = 0
    tcnt = 0
    for Ridx in range(NUM_FLX):
        R = selected_reactions[Ridx]
        OUTPUT_NAMES.append(R)
        MIN_STAT = MAX_STAT = 0
        # MIN
        ##  cplx_setObjective(cpx, R, expr=None, sense='min', reset=True)
        # TODO: bgoli: see whether this also works with 'minimize'
        cplx_setObjective(cpx, 'min%s' % R, [(1, R)], 'min', reset=True)
        ##  cplx_setBounds(c, id, min=None, max=None) # think about this
        MIN_STAT = cplx_Solve(cpx, method='o')
        if MIN_STAT == 1: # solved
            min_sol, min_oid, min_oval = cplx_getOptimalSolution(cpx)
        elif MIN_STAT == 2: # unbound
            min_sol = None
            min_oid = None
            min_oval = -numpy.Inf            
        else: # other failure
            min_sol = None
            min_oid = None
            min_oval = numpy.NaN
        if debug:
            cplx_writeLPtoLPTfile(cpx, '%smin' % R, title='min%s=%s' % (R,min_oval), Dir=debug_dir)

        # MAX
        ##  cplx_setObjective(cpx, R, expr=None, sense='max', reset=True)
        cplx_setObjective(cpx, 'max%s' % R, expr=None, sense='max', reset=False)
        ##  cplx_setBounds(c, id, min=None, max=None) # think about this
        MAX_STAT = cplx_Solve(cpx, method='o')
        if MAX_STAT == 1: # solved
            max_sol, max_oid, max_oval = cplx_getOptimalSolution(cpx)
        elif MAX_STAT == 2: # unbound
            max_sol = None
            max_oid = None
            max_oval = numpy.Inf
        else: # other fail
            max_sol = None
            max_oid = None
            max_oval = numpy.NaN            
        if debug:
            cplx_writeLPtoLPTfile(cpx, '%smax' % R, title='max%s=%s' % (R,max_oval), Dir=debug_dir)

        # print 'DEBUGmin:', type(min_sol), min_oid, min_oval, MIN_STAT
        # print 'DEBUGmax:', type(max_sol), max_oid, max_oval, MAX_STAT

        OUTPUT_ARRAY[Ridx,0] = pre_sol[R]
        if REDUCED_COSTS.has_key(R):
            OUTPUT_ARRAY[Ridx,1] = REDUCED_COSTS[R]
        OUTPUT_ARRAY[Ridx,2] = min_oval
        OUTPUT_ARRAY[Ridx,3] = max_oval
        OUTPUT_ARRAY[Ridx,4] = abs(max_oval - min_oval)
        OUTPUT_ARRAY[Ridx,5] = MIN_STAT
        OUTPUT_ARRAY[Ridx,6] = MAX_STAT
        if MAX_STAT != 1 or MIN_STAT != 1:
            print 'Solver fail for reaction', R
            print OUTPUT_ARRAY[Ridx]
            print 'Minimization status:', MIN_STAT
            print 'Maximization status:', MAX_STAT
            ##  raw_input(R)
            time.sleep(0.1)
        ##  if R == 'R_NDPK1':
            ##  raw_input('R_NDPK1')        
        ##  if R == 'R_ACt2rpp':
            ##  raw_input('R_ACt2rpp')
        cntr += 1
        if cntr == 200:
            tcnt += cntr
            print 'FVA has processed %s of %s reactions' % (tcnt, NUM_FLX)
            cntr = 0
    if quiet:
        cplx_setOutputStreams(cpx, mode='default')

    print 'Output array has columns:'
    print 'Reaction, Reduced Costs, Variability Min, Variability Max, abs(Max-Min), MinStatus, MaxStatus'
    return OUTPUT_ARRAY, OUTPUT_NAMES



def cplx_MultiFluxVariabilityAnalysis(lp, selected_reactions=None, tol=None, rhs_sense='lower', optPercentage=100.0, work_dir=None, debug=False):
    """
    Perform a flux variability analysis on a multistate LP

     - *lp* a multistate LP
     - *selected reactions* [default=None] means use all reactions otherwise use the reactions listed here
     - *pre_opt* [default=True] attempt to presolve the FBA and report its results in the ouput
     - *tol*  [default=None] do not floor/ceiling the objective function constraint, otherwise round of to *tol*
     - *rhs_sense* [default='lower'] means objC >= objVal the inequality to use for the objective constraint can also be *upper* or *equal*
     - *optPercentage* [default=100.0] means the percentage optimal value to use for the RHS of the objective constraint: optimal_value*(optPercentage/100.0)
     - *work_dir* [default=None] the FVA working directory for temporary files default = cwd+fva
     - *debug* [default=False] if True write out all the intermediate FVA LP's into work_dir

    and returns an array with columns::

     Reaction, Reduced Costs, Variability Min, Variability Max, abs(Max-Min), MinStatus, MaxStatus

    and a list containing the row names.
    """
    assert RuntimeError, '\nThis method only works on the L1 distance and has been removed for now'
    
    gc.enable()
    pre_opt = True
    debug_dir = work_dir
    if debug_dir == None:
        debug_dir = os.path.join(os.getcwd(), 'fva')
    else:
        assert os.path.exists(debug_dir), '\nWhat did you think would happen now!'
        debug_dir = os.path.join(debug_dir, 'fva')
    if debug and not os.path.exists(debug_dir):
        os.mkdir(debug_dir)

    OPTIMAL_PRESOLUTION = None
    pre_sol = pre_oid = pre_oval = None
    REDUCED_COSTS = {}
    if pre_opt:
        cplx_Solve(lp, method='o')
        if cplx_getSolutionStatus(lp) == 'LPS_OPT':
            print 'Valid FVA presolution'
            OPTIMAL_PRESOLUTION = True
            pre_sol, pre_oid, pre_oval = cplx_getOptimalSolution(lp)
            REDUCED_COSTS = cplx_getReducedCosts(lp)
        else:
            print 'Invalid FVA presolution'
            OPTIMAL_PRESOLUTION = False
            pre_sol = {}
            pre_oval = 0.0
            pre_oid = 'DisMinObj'
            raise RuntimeError, '\nFVA presolve failed to optimize this model and cannot continue!'
    ##  else:
        ##  pre_sol = {}
        ##  for r in fba.reactions:
            ##  pre_sol.update({r.getPid : 0.0})
        ##  pre_oval = fba.objectives[fba.activeObjIdx].value
        ##  pre_oid = fba.objectives[fba.activeObjIdx].getPid()
    ##  for r in fba.reactions:
        ##  r.value = pre_sol[r.id]

    # generate new constraint from old objective value (use non-zero coefficients)
    LCS = lp.objective.get_linear()
    LCN = lp.variables.get_names()
    new_constraint = []
    for c in range(len(LCS)):
        if LCS[c] != 0.0:
            new_constraint.append((LCS[c], LCN[c]))

    del LCS, LCN

    # set objective constraint
    if rhs_sense == 'equal':
        cplx_setSingleConstraint(lp, 'ObjCstr', expr=new_constraint, sense='E', rhs=pre_oval)
        ##  cplx_setSingleConstraint(lp, 'ObjCstr', expr=[(1, pre_oid)], sense='E', rhs=pre_oval)
    elif rhs_sense == 'upper':
        if tol != None:
            ub =  numpy.ceil(pre_oval/tol)*tol*optPercentage/100.0
        else:
            ub = pre_oval*(optPercentage/100.0)
        cplx_setSingleConstraint(lp, 'ObjCstr', expr=new_constraint, sense='L', rhs=ub)
        ##  cplx_setSingleConstraint(lp, 'ObjCstr', expr=[(1, pre_oid)], sense='L', rhs=ub)
    elif rhs_sense == 'lower':
        if tol != None:
            lb =  numpy.floor(pre_oval/tol)*tol*optPercentage/100.0
        else:
            lb = pre_oval*(optPercentage/100.0)
        cplx_setSingleConstraint(lp, 'ObjCstr', expr=new_constraint, sense='G', rhs=lb)
        ##  cplx_setSingleConstraint(lp, 'ObjCstr', expr=[(1, pre_oid)], sense='G', rhs=lb)

    if debug:
        cplx_writeLPtoLPTfile(lp, 'fva_base_%s' % time.time(), title=None, Dir=debug_dir)

    NUM_FLX = lp.variables.get_num()
    VARIABLE_NAMES = lp.variables.get_names()
    print 'Number of variables: %s' % NUM_FLX
    if selected_reactions != None:
        for r in selected_reactions:
            assert r in VARIABLE_NAMES, "\n%s is not a valid reaction name" % r
        NUM_FLX = len(selected_reactions)
        print 'Number of user selected variables: %s' % NUM_FLX
    try:
        OUTPUT_ARRAY = numpy.zeros((NUM_FLX, 7), numpy.double)
        ##  OUTPUT_ARRAY = numpy.zeros((NUM_FLX, 7))
    except AttributeError:
        OUTPUT_ARRAY = numpy.zeros((NUM_FLX, 7))
    OUTPUT_NAMES = []
    iter_cnt = 0
    for Ridx in range(NUM_FLX):
        if selected_reactions != None:
            R = selected_reactions[Ridx]
        else:
            R = VARIABLE_NAMES[Ridx]
        iter_cnt += 1
        OUTPUT_NAMES.append(R)
        MIN_STAT = MAX_STAT = 0
        # MIN
        ##  cplx_setObjective(lp, R, expr=None, sense='min', reset=True)
        cplx_setObjective2(lp, 'min%s' % R, VARIABLE_NAMES, [(1, R)], 'min', reset=True)
        ##  cplx_setBounds(c, id, min=None, max=None) # think about this
        cplx_Solve(lp, method='o')
        if cplx_getSolutionStatus(lp) == 'LPS_OPT':
            #DBG# min_sol, min_oid, min_oval = cplx_getOptimalSolution(lp)
            min_sol, min_oid, min_oval = cplx_getOptimalSolution2(lp, VARIABLE_NAMES)
            ##  min_sol, min_oid, min_oval = 0.0,0.0,0.0
            MIN_STAT = 1
        else:
            ##  min_sol, min_oid, min_oval = 0.0
            min_sol = None
            min_oid = None
            min_oval = numpy.NaN
        if debug:
            cplx_writeLPtoLPTfile(lp, '%smin' % R, title='min%s=%s' % (R,min_oval), Dir=debug_dir)

        # MAX
        ##  cplx_setObjective(lp, 'max%s' % R, expr=None, sense='max', reset=False)
        cplx_setObjective2(lp, 'max%s' % R, VARIABLE_NAMES, expr=None, sense='max', reset=False)
        ##  cplx_setBounds(c, id, min=None, max=None) # think about this
        cplx_Solve(lp, method='o')
        if cplx_getSolutionStatus(lp) == 'LPS_OPT':
            #DBG# max_sol, max_oid, max_oval = cplx_getOptimalSolution(lp)
            max_sol, max_oid, max_oval = cplx_getOptimalSolution2(lp, VARIABLE_NAMES)
            ##  max_sol, max_oid, max_oval = 0.0,0.0,0.0
            MAX_STAT = 1
        else:
            ##  max_sol, max_oid, max_oval = 0.0
            max_sol = None
            max_oid = None
            max_oval = numpy.NaN
        if debug:
            cplx_writeLPtoLPTfile(lp, '%smax' % R, title='max%s=%s' % (R,max_oval), Dir=debug_dir)

        OUTPUT_ARRAY[Ridx,0] = pre_sol[R]
        if REDUCED_COSTS.has_key(R):
            OUTPUT_ARRAY[Ridx,1] = REDUCED_COSTS[R]
        OUTPUT_ARRAY[Ridx,2] = min_oval
        OUTPUT_ARRAY[Ridx,3] = max_oval
        OUTPUT_ARRAY[Ridx,4] = abs(max_oval - min_oval)
        OUTPUT_ARRAY[Ridx,5] = MIN_STAT
        OUTPUT_ARRAY[Ridx,6] = MAX_STAT
        del max_sol, max_oid, max_oval, min_sol, min_oid, min_oval
        if MIN_STAT == 0.0 or MIN_STAT == 0.0:
            print 'Solver fail for reaction'
            print R, OUTPUT_ARRAY[Ridx]
        if iter_cnt == 100:
            iter_cnt = 0
            ##  mdtempFile = os.path.join(debug_dir,'_memDump.lp')
            ##  lp.write(mdtempFile, filetype='lp')
            ##  del lp
            gc.collect()
            ##  reload(cplex)
            print 'Collecting garbage'
            ##  time.sleep(2)
            ##  lp = cplex.Cplex(mdtempFile)
    for r in OUTPUT_ARRAY:
        print r
    print 'Output array has columns:'
    print 'Reaction, Reduced Costs, Variability Min, Variability Max, abs(Max-Min), MinStatus, MaxStatus'
    del lp
    return OUTPUT_ARRAY, OUTPUT_NAMES


def cplx_WriteFVAtoCSV(id, fva, names, Dir=None, fbaObj=None):
    """
    Takes the resuls of a FluxVariabilityAnalysis method and writes it to a nice
    csv file. Note this method has been refactored to `CBWrite.WriteFVAtoCSV()`.
    
     - *id* filename_base for the CSV output
     - *fva* FluxVariabilityAnalysis() OUTPUT_ARRAY
     - *names* FluxVariabilityAnalysis() OUTPUT_NAMES
     - *Dir* [default=None] if set the output directory for the csv files
     - *fbaObj* [default=None] if supplied adds extra model information into the output tables
     
    """
    
    print '\nThis method is now a wrapper for the newer CBWrite.WriteFVAtoCSV method, provided for back compatbility.\n'
    
    CBWrite.WriteFVAtoCSV(id, fva, names, Dir, fbaObj)
